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60.1511 



DETERMINATION OF STRESS CHARACTERISTICS 
OF EARTH FORMATIONS 

FIELD OF THE INVENTION 

This invention relates to investigation of earth formations and, more 
particularly, to a method and apparatus for obtaining properties of earth 
formations using sonic logging and determining stress characteristics of the 
earth formations. 

BACKGROUND OF THE INVENTION 

It is well known that mechanical disturbances can be used to establish 
acoustic waves in earth formations surrounding a borehole, and the properties 
of these waves can be measured to obtain important information about the 
formations through which the waves have propagated. Parameters of 
compressional, shear, and Stoneley waves, such as their velocity (or its 
reciprocal, slowness) in the formation and in the borehole, can be indicators of 
formation characteristics that help in evaluation of the location and/or 
producibility of hydrocarbon resources. 

An example of a logging device that has been used to obtain and 
analyze sonic logging measurements of formations surrounding an earth 
borehole is called a Dipole Shear Sonic Imager ("DSI" - trademark of 
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Schlumberger), and is of the general type described in Harrison et al., 
"Acquisition and Analysis of Sonic Waveforms From a Borehole Monopole 
And Dipole Source For The Determination Of Compressional And Shear 
Speeds And Their Relation To Rock Mechanical Properties And Surface 
Seismic Data", Society of Petroleum Engineers, SPE 20557, 1990. 

An acoustic source in a fluid-filled borehole generates 
headwaves as well as relatively stronger borehole-guided modes. 
A standard sonic measurement system consists of placing a 
piezoelectric or electrodynamic source and an array of 
hydrophone receivers inside a fluid-filled borehole. The acoustic 
source is configured in the form of either a monopole or a dipole 
source. The source bandwidth typically ranges from 0.5 to 20 
kHz. A monopole source generates primarily the lowest-order 
axisymmetric mode, also referred to as the Stoneley mode, 
together with compressional and shear headwaves. In contrast, a 
dipole source primarily excites the lowest-order borehole flexural 
mode together with compressional and shear headwaves. The 
headwaves are caused by the coupling of the transmitter acoustic 
energy to plane waves in the formation that propagate along the 
borehole axis. An incident compressional wave in the borehole 
fluid produces critically refracted compressional waves in the 
formation. These refracted waves along the borehole surface are 
known as compressional headwaves. The critical incidence angle 
0i = sin" 1 (Vf/V c ), where V f is the compressional wave speed in the 
borehole fluid; and V c is the compressional wave speed in the 
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formation. As the compressional headwave travels along the 
interface, it radiates energy back into the fluid that can be 
detected by hydrophone receivers placed in the fluid-filled 
borehole. In fast formations, the shear headwave can be 
similarly excited by a compressional wave at the critical incidence 
angle 0i = sin" 1 (V f /V c ), where V f is the shear wave speed in the 
formation. In a homogeneous and isotropic model of fast 
formations, compressional and shear headwaves can be 
generated by a monopole source placed in a fluid-filled borehole 
for determining the formation compressional and shear wave 
speed (Kimball and Marzetta, 1984). It is known that refracted 
shear headwaves cannot be detected in slow formations (where 
the shear wave velocity is less than the borehole-fluid 
compressional velocity) with receivers placed in the borehole 
fluid. In slow formations, formation shear velocities are obtained 
from the low-frequency asymptote of borehole flexural dispersion. 
Standard processing techniques exist for the estimation of 
formation shear velocities in either fast or slow formations from 
an array recorded dipole waveforms. The radial depth of 
investigation in the case of headwave logging is dependent on the 
transmitter to receiver spacing and any in-situ radial variations in 
formation properties that might be present. The radial depth of 
investigation in the case of modal logging is well characterized 
and it extends to about a wavelength. This implies that the low- 
and high-frequencies probe deep and shallow into the formation, 
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respectively. Existing sonic tools, such as the above-referenced 
"DSI" are capable of providing the far-field compressional, fast- 
shear, and slow-shear velocities away from the near-wellbore 
alterations. 

The effective formation volume probed by borehole sonic 
waves is generally, subject to (1) near-wellbore stress 
distributions caused by the far-field formation stresses; (2) 
wellbore mud pressure; and (3) pore pressures. Formation 
stresses play an important role in geophysical prospecting and 
developments of oil and gas reservoirs. Both the direction and 
magnitude of these stresses are required in (a) planning for 
borehole stability during directional drilling, (b) hydraulic 
fracturing for enhanced production, and (c) selective perforation 
for prevention of sanding during production. A detailed 
knowledge of formation stresses also helps in management of 
reservoirs that are prone to subsidence caused by a significant 
reduction in pore pressure and an associated increase in the 
effective stress that exeeds the in-situ rock strength. In addition, 
the magnitude and orientation of the in-situ stresses in a given 
field have a significant influence on permeability distribution that 
help in a proper planning of wellbore trajectories and injection 
schemes for water or steam flooding. As estimates of stresses 
from borehole measurements are improved, it is not uncommon to 
find that the regional tectonic stress model involving large global 
averages is significantly different than the local stresses around a 
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borehole that affects the reservior producibility and near-wellbore 
stability. 

The formation stress state is characterized by the 
magnitude and direction of the three principal stresses. 
Generally, the overburden stress yields the principal stress in the 
vertical direction. The overburden stress (S v ) is reliably obtained 
by integrating the formation mass density from the surface to the 
dept of interest. Consequently, estimating the other two principal 
stresses (S Hm ax and S hm in) in the horizontal plane is the remaining 
task necessary to fully characterize the formation stress state. 
Figure 1 shows various applications of the formation principal 
stresses together with the well bore (P p ) pressures in well 
planning, well bore stability and reservoir management. 

The far-field formation stresses can be expressed in terms 
of three principal stresses - they are referred to as the 
overburden stress S v , maximum horizontal stress SHmax, and 
minimum horizontal stress S hm in as shown in Figure 2. It is known 
that the wellbore mud overpressure has rather small effect on 
borehole flexural velocities at low frequencies. In contrast, 
flexural velocities at low frequencies are more sensitive to 
changes in the far-field formation horizontal stresses. 

Sonic velocities in formations change as a function of rock 
lithology/mineralogy, porosity, clay content, fluid saturation, 
stresses, and temperature. To estimate changes in the formation 
stress magnitudes from measured changes in sonic velocities, it 
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is necessary to select a depth interval with a reasonably uniform 
lithology, clay content, saturation, and temperature so that the 
measured changes in velocities can be related to corresponding 
changes in formation stress magnitudes. Any change in porosity 
in the chosen depth interval is accounted for by corresponding 
changes in the formation effective bulk density and stiffnesses. 
Assuming that the measured changes in sonic velocities are 
largely caused by changes in stress magnitudes, it is possible to 
invert borehole sonic velocities for the estimation of changes in 
formation stress magnitudes (see U.S. Patent No. 6,351,991). 
These techniques neither require the presence of a wellbore 
failure in the form of a breakout or fractures, nor the in-situ rock 
tensile and compressive strengths. The technique is based on 
changes in velocities of multiple sonic waves caused by changes 
in prestress in the propagating medium. It is important to 
ascertain that one uses sonic velocities outside any mechanically 
altered annulus caused by borehole-induced stress 
concentrations and plastic deformations that might also occur. 

The propagation of small amplitude waves in homogeneous 
and anitsotropic materials is described by the linear equations of 
motion. However, when the material is prestressed, the 
propagation of such waves are properly described by equations of 
motion for small dynamic fields superposed on a static bias. A 
static bias represents any statically deformed state of the medium 
caused by an externally applied load or residual stresses. The 
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equations of nonlinear elasticity (see B.K. Sinha, "Elastic Waves 
In Crystals Under A Bias", Ferroelectrics, (41), 61-73, 1982; A.N. 
Norris, B.K. Sinha, and S. Kostek, "Acoustoelasticity Of 
Solid/Fluid Composite Systems", Geophys. J. Int., (118), 439-446, 
1994). The linear equations of motion for isotropic materials 
contain two independent stiffnesses that are defined by the 
dynamic Young's modulus (Y) and Poisson's ratio (v) or 
equivalently, by the two Lame parameters, (A and p). The 
equations of motion or prestressed isotropic materials contain 
three additional nonlinear stiffness constants (Cm, c 14 4, and 
C155). 
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SUMMARY OF THE INVENTION 



Embodiments of the present invention include new 
techniques for the estimation of formation horizontal stresses. 
One of the techniques includes a multi-frequency inversion of 
cross-dipole dispersions for both the maximum and minimum 
horizontal stress magnitudes. A second technique involves 
inverting flexural velocity-dispersion gradients for estimating the 
isotropic horizontal stress in depth intervals with a reasonably 
uniform lithology. The second technique is, generally, applicable 
to shale intervals where dipole dispersion crossovers are not 
observed, implying isotropic horizontal stresses. 

In accordance with an embodiment of the invention, a 
method is set forth for determining the maximum and minimum 
horizontal stresses of formations surrounding a borehole, 
including the following steps: suspending a logging device in the 
borehole; transmitting sonic energy from the logging device to 
establish flexural waves in the formations; receiving, at the 
logging device, sonic energy from the flexural waves, and 
producing, from the received sonic energy, measurement signals 
at a number of frequencies; determining, at said number of 
frequencies, the fast and slow flexural wave velocities in the 
formations, to obtain fast and slow flexural velocity dispersions; 
establishing a model of formation stresses in which stresses of a 
loaded state are represented by the sum of an omnidirectional 
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hydrostatically loaded mean reference stress, a vertical stress 
perturbation, and maximum and minimum horizontal stress 
perturbations; establishing an inversion model that includes 
inputs from the fast and slow flexural velocity dispersions and 
also includes unknown horizontal perturbations of said model of 
formation stresses; deriving, from the inversion model, estimates 
of the maximum and minimum horizontal stress perturbations; and 
determining, from the estimates of the maximum and minimum 
horizontal stress perturbations and the mean reference stress, 
estimates of the maximum and minimum horizontal stresses of the 
formations. 

A preferred form of this embodiment of the invention further 
includes the step of determining stress coefficients of the fast 
and slow flexural velocities as a function of frequency, and using 
the determined stress coefficients to obtain an estimate of said 
vertical stress perturbation. This form of the invention further 
includes deriving a reference flexural velocity dispersion, and the 
determination of the estimate of the vertical stress perturbation 
also depends on said reference flexural velocity dispersion. In 
this form, the inversion model includes differences between the 
fast and slow flexural velocity dispersions, normalized by the 
reference flexural velocity dispersion. The reference flexural 
velocity dispersion is derived from an assumed isotropic and 
homogeneous formation model which is, in turn, formulated using 
using measured monopole compresssional and shear velocities. 
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In accordance with a further embodiment of the invention, a 
method is set forth for determining isotropic horizontal stresses in 
a shale interval of formations surrounding a borehole, including 
the following steps: (a) suspending a logging device in the 
borehole; (b) transmitting sonic energy from the logging device to 
establish flexural waves in the formations; (c) receiving, at the 
logging device, sonic energy from the flexural waves, and 
producing, from the received sonic energy, measurement signals 
at a number of frequencies; (d) determining, at said number of 
frequencies, the flexural wave velocities in the formations, to 
obtain a flexural velocity dispersion; (e) performing the steps (a) 
through (d) at vertically spaced apart upper and lower depth 
levels in the shale region to obtain upper and lower flexural 
velocity dispersions; (f) establishing models of formation stresses 
at the upper and lower depth levels in which stresses of a loaded 
state are represented by the sum of an omnidirectional 
hydrostatically loaded mean reference stress, a vertical stress 
perturbation, and a horizontal stress perturbation; (g) establishing 
inversion models that include inputs from the upper and lower 
flexural velocity dispersions and also includes unknown vertical 
and horizontal perturbations of the model of formation stresses at 
the upper and lower depth levels; (h) deriving, from the inversion 
models, estimates of the vertical and horizontal stress 
perturbations at the upper and lower depth levels; and (i) 
determining, from the estimates of the vertical and horizontal 
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stress perturbations at the upper and lower depth levels and the 
mean reference stress, estimates of the horizontal stress at the 
upper and lower depth levels in the shale interval of the 
formations. 

Further features and advantages of the invention will 
become more readily apparent from the following detailed 
description when taken in conjunction with the accompanying 
drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 
Figure 1 is a diagram illustrating applications of formation 

stresses in well planning, wellbore stability and reservoir 

management. 

Figure 2 is a schematic diagram of a borehole in the 

presence of formation principal stresses. The three principal 

stresses are assumed to be different. 

Figure 3 is a diagram, partially in block form, of a type of 

apparatus that can be used in practicing embodiments of the 

invention. 

Figure 4 is a simplified diagram of a type of downhole 
logging device that can be used in practicing embodiments of the 
invention. 

Figure 5 is a diagram illustrating, in cross-section, the 
placement of hydrophones that can be used at a receiver station 
in the logging device of Figure 4. 

Figure 6 is a block diagram of a portion of the electronics of 
the Figure 4 logging device. 

Figure 7 illustrates a model of decomposition of formation 
effective stresses in the far-field away from the borehole into a 
hydrostatically loaded reference and perturbations of triaxial 
stresses. 

Figure 8 is a graph of an example of the fast-, slow-, and 
reference flexural velocity dispersions. Normalized differences 
between the fast- and slow-flexural velocities at various 
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wavenumbers are input to the inversion model for the estimation 
of formation nonlinear constants and Ao H and Aa h . kj is the axial 
wavenumber, and a is the borehole radius. The index i denotes 
the number of input velocity Vi and frequency fj used in the 
inversion model. 

Figure 9 is a graph of an example of the fast-, slow-, and 
reference flexural velocity dispersions. Normalized differences 
between (1) the fast-flexural and reference flexural velocities; 
and (2) the slow-flexural and reference flexural velocities at 
various wavenumbers are input to the inversion model for the 
estimation of Aov- The notation is the same as in Figure 8. 

Figure 10 is a schematic diagram of a borehole in the 
presence of formation principal stresses. The horizontal stresses 
are assumed to the same in this case. 

Figure 11 illustrates a model of decomposition of formation 
effective streses in the far-field away from the borehole into a 
hydrostatically loaded reference and perturbations of triaxial 
stresses. In this case, stresses in the horizontal plane 
are the same. 

Figure 12 is a graph of an example of the bottom-, top-, and 
reference flexural velocity dispersions. Normalized differences 
between the bottom- and top-flexural velocities at various 
wavenumbers are input to the inversion model for the estimation 
of formation nonlinear constants and ( Act* -Act*). The notation is 
the same as in Figure 8. 
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Figure 13 is a graph of an example of the bottom-, top-, and 
reference flexural velocity dispersions. Normalized differences 
between (1) the bottom-flexural and reference flexural velocities; 
and (2) the top-flexural and reference flexural velocities at 
various wavenumbers are input to the inversion model for the 
estimation of Ao H . The notation is the same as in Figure 8. 

Figures 14A and 14B, and 14C, when placed one beneath 
another, show a flow diagram of a routine for programming a 
processor to implement an embodiment of the invention. 

Figures 15A and 15B, and 15C, when placed one beneath 
another; show a flow diagram of a routine for programming a 
processor to implement another embodiment of the invention. 
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DETAILED DESCRIPTION 



Referring to Figure 3, there is shown a type of apparatus which can be 
used in practicing embodiments of the invention. Subsurface formations 231 
are traversed by a borehole 232 which is typically, although not necessarily, 
filled with drilling fluid or mud. A logging tool 210 is suspended on an armored 
cable 212 and may have optional centralizers (not shown). The cable 212 
extends up the borehole, over a sheave wheel 220 on a derrick 221 to a 
winch forming part of surface equipment 250. Known depth gauging 
apparatus (not shown) is provided to measure cable displacement over the 
sheave wheel 220 and accordingly the depth of the logging tool 210 in the 
borehole 232. A device of a type well known in the art is included in the tool 
210 to produce a signal indicative of orientation of the body of the tool 210. 
Processing and interface circuitry within the tool 210 amplifies, samples and 
digitizes the tool's information signals for transmission and communicates 
them to the surface equipment 250 via the cable 212. Electrical power and 
control signals for coordinating operation of the tool 210 are generated by the 
surface equipment 250 and communicated via the cable 212 to circuitry 
provided within the tool 210. The surface equipment includes processor 
subsystem 270 (which can typically include a microprocessor, memory, clock 
and timing, and input/output functions - not separately shown), standard 
peripheral equipment (not separately shown), and recorder 226. 

The logging device 210 may be, for example, of a type first noted 
above and known as a Dipole Shear Sonic Imager ("DSI" - trademark of 
Schlumberger), generally described in Harrison et al., "Acquisition and 
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Analysis of Sonic Waveforms From a Borehole Monopole and Dipole Source 
for the Determination of Compressional and Shear Speeds and Their Relation 
to Rock Mechanical Properties and Surface Seismic Data", Society of 
Petroleum Engineers, SPE 20557, 1990. It will be understood, however, that 
any suitable logging device can be utilized. Further details of the logging 
device 210 of this example are shown in Figure 4. The logging device 210 
includes crossed dipole transmitters 315 and 320 (only one end of dipole 320 
being visible) and a monopole transmitter 325, so that waves including 
compressional, shear, Stoneley, and flexural can be excited. Orthogonal 
transmitter dipoles are preferred in an embodiment hereof that employs two 
pairs of transmitter dipoles, but transmitter dipoles at a known non-zero angle 
that is not a right angle can also be used. As seen in the illustration, the DSI 
orthogonal dipole transmitters of this example are not collocated, although 
they could be. The present invention can make use of either non-collocated 
or collocated orthogonal dipole transmitters without any loss of generality. 
Eight, or other suitable number, of spaced apart receiver stations, designated 
331 through 338 each comprise four receiver hydrophones mounted 
azimuthally at ninety degree intervals in the surface of the cylindrical logging 
device. Figure 5 shows the hydrophones, designated A, B, C, and D. In an 
example shown in Figure 6, an X component can be obtained by subtracting 
the signals received at A and C (i.e., A-C), and a Y component can be 
obtained by subtracting the signals received at B and D (i.e., B-D). With four 
receiver elements at each receiver station, there are a total of thirty-two 
receiver elements in this example. The receiver stations are also configurable 
for monopole reception. 
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The transmitter electronics contain a power amplifier and switching 
circuitry capable of driving the two crossed-dipole transmitter elements and 
the monopole element from a programmable waveform. Separate waveforms 
with appropriate shape and frequency content can be used for dipole, 
Stoneley and compressional measurements. The receiver electronics 
processes the signals from the 32 individual receiver elements located at the 
eight receiver stations which are spaced six inches apart. At each station, 
four receivers are mounted as shown in Figure 5 which allows measurement 
of the dipole and crossed-dipole waveforms by differencing the outputs from 
opposite receivers, as previously described. Summing the outputs of the 
receivers can be used to produce a monopole equivalent signal. As further 
described in Harrison et al., supra, the receiver electronics multiplexers, 
filters, amplifies and channels the signals from the 32 receiver elements to 
eight parallel signal paths. These eight parallel analog signals are passed to 
an acquisition electronics cartridge where eight 12-bit analog-to-digital 
converters digitize the signals from the receiver electronics. The telemetry 
circuitry passes the digitized information to the earth's surface. 

Figure 6 shows an example of the acquisition signal path in block 
diagram for for one of the eight (or other suitable number of) receiver stations, 
as described in Harrison et al., supra. Each receiver has its own charge 
preamplifier (represented at 505). The output of the receivers, odd numbered 
pairs being in-line with the upper dipole transmitter and even numbered pairs 
with the lower dipole transmitter, passes into both a summing circuit (for 
monopole measurements) and a differencing circuit (for dipole 
measurements), as represented at 510. Under software control the sum or 
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difference is selected by a multiplexer stage (block 520) and the signal 
passed to one of eight programmable gain amplifier stages (540) and filters 
(545). The other similar channels are represented by block 550. The eight 
parallel analog signals are passed to eight parallel 12-bit A/D converters 
(represented at 560) where simultaneous waveform digitization is performed. 
If desired, more bits can, of course, be used to advantage. After digitization, 
the eight waveforms are passed to the memory section associated with 
downhole processor 580. The processor also provides control signals and 
waveforms to transmitter and receiver electronics. An alternate path directs 
the eight analog receiver signals into threshold crossing detection circuitry or 
digital first motion detection, as represented at block 565. This circuitry 
detects the time of all up or down going threshold crossings. The digitized 
waveform data and the threshold crossing time data are passed to the surface 
using telemetry circuitry 590. It will be understood that more advanced tool 
implementations, having additional transmitters, receivers, and/or transmitter- 
to-receiver (T/R) spacings, and more powerful processing capabilities, can be 
used even more advantageously, consistent with the principles hereof, in 
practicing embodiments of the invention. 

In the Figure 4 embodiment, the processing of signals recorded uphole 
can be implemented using a processor 270, such as a suitably programmed 
general purpose digital processor with memory and peripherals conventionally 
provided. It will be understood, however, that the processing need not be 
performed at the wellsite, and that signals derived at the wellsite can be 
processed at a remote location. It will also be understood that other suitable 
logging tools can be employed in practicing the invention. 
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Elastic wave velocities in porous materials change as a function of 
effective stresses in the propagating medium. These effective stresses a ti are 
defined (see Nur and Byerlee, "An Exact Effective Stress Law For Elastic 
Deformation of Rock With Fluids", Journal of Geophysical Research, vol. 76, 
pp. 6414-6419, 1971) by 

<r 9 -S,-a8,P,, (1) 
where Sy is the Kronecker delta and the Biot parameter a is given by 

a=1 -f' ( 2 ) 

where K is the bulk modulus of the dry aggregate and K s is the intrinsic bulk 
modulus of the solid matrix. Even though the porosity effect is not explicit in 
this expression, it is included in the value of the effective bulk modulus K of 
the dry aggregate. In the case of a nonporous and impermeable formation, K 
« K s , and a = 0. In contrast, in natural soils, K « Ks, and a = 1 . In porous 
rocks, the Biot parameter is affected by both the porosity, pore shape, and 
connectivity (see Walsh, 'The Effects of Cracks on The Compressibility of 
Rocks", Journal of Geophysical Research, vol. 70, p. 381 , 1965). Generally, 
<J> < a < 1 , in porous rocks with porosity <t>. More importantly, the Biot 
parameter a n in a normally compacted shale formation with open pores is 
typically more than that in an overpressured (a G ) shale with larger effective 
porosity and closed pores, i.e., a n > a Q . 

Consider a vertical fluid-filled borehole with its axis parallel to the Xr 
direction and the maximum and minimum horizontal stresses parallel to the 
X 2 - and X 3 -directions, respectively. Triaxial formation stresses with a vertical 
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overburden stress as one of the principal stresses can be decomposed into a 
hydrostatically loaded isotropic reference and perturbations in the three 
principal stresses Ao v , Ao H , and Ao h , as shown in Figure 7. Note that the 
mean stress o vo in the isotropic reference state is not known at this point. 
However, this assumed state is defined from the measured compressional 
and monopole shear velocities at the chosen depth so that small perturbations 
in the three principal stresses Aa v , Aa H , and Ao h , would lead to the actual in- 
situ stresses at this depth. All of these stresses in Figure 7 are far-field 
stresses beyond any stress-concentration annulus caused by the presence of 
a borehole. 

When the propagating medium is prestressed, the propagation of 
small-amplitude waves are properly described by equations of motion for 
small dynamic fields superposed on a prestressed state. A prestressed state 
represents any statically deformed state of the medium caused by an 
externally applied load or residual stresses. The equations of motion for small 
dynamic fields superposed on a static bias are derived from the rotationally 
invariant equations of nonlinear elasticity by making a Taylor expansion of the 
quantities for the dynamic state about their values in the biasing (or 
intermediate) state (see N. Norris, B.K. Sinha, and S. Kostek, 
"Acoustoelasticity Of Solid/Fluid Composite Systems, Geophys. J. Int., (118), 
439-446, 1994; B.K. Sinha and S. Kostek, "Stress-Induced Azimuthal 
Anisotropy In Borehole Flexural Waves", Geophysics, (61) 6, 1899-1907, 
1996). 

When the biasing state is inhomogeneous, the effective elastic 
constants are position dependent, and a direct solution of the boundary value 
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problem is not possible. In this situation, a perturbation procedure can readily 
treat spatially varying biasing states, such as those in the case of nonuniform 
radial and circumferential stress distributions away from the borehole, and the 
corresponding changes in the Stoneley and flexural-wave velocities can be 
calculated as a function of frequency (Norris et al., 1994, supra; Sinha et al., 
1 995, supra). Referenced to the statically deformed state of the formation (or 
the intermediate configuration), the modified Piola-Kirchhoff stress tensor P ai 

as defined in Norris et al. 1994, supra, can be employed in a perturbation 
model to obtain the following expression for the first-order perturbation in the 
eigenfrequency oj m for a given wavenumber k z , 



Acq _ v 



°>» 2aj 2 m j Po uyyv 



(3) 



where 



AP^=H^u^, (4) 
Hqjyp = >W +Vy r + **** << S A - S «t S jA ( 5 ) 

hcqyp ~ ~ C cgyp W S,S + C cg^AB^AB 

+ w r,r c *i Pfi + w fi,Q c cM> (6) 

T afi = c ap r s w sr. ( 7 ) 

E AB =^(w AB +w BA ). (8) 

Equations (3)-(8) use the Cartesian tensor notion and the convention that a 
comma followed by an index P denotes differentiation with respect to the 
spatial coordinate X p . The summation convention for repeated tensor indices 
is also implied. Although h ajY p exhibits the usual symmetries of the second- 
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order constants of linear elasticity, the effective elastic stiffness tensor H a j Y p 
does not have these properties as is evident from equation (5). Before 
discussing the various quantities in equations (3)-(8), it is noted that the 
present position of material points may be written as 

y(X,t) = X + w(X) + u(X,t), (9) 
where w denotes the displacement caused by the applied static loading of 
material points with position vector X in the assumed reference isotropic state 
(also defined as the reference state), and u denotes the small dynamic 
displacement vector of material points above and beyond that due to the static 
deformation. The small field Piola-Kirchhoff stress P aj in the intermediate 
state can be decomposed into two parts 

P qi =P^+AP^ (10) 

where 

(11) 

with AP aj being defined by equations (4)-(8), and the superscript "L" denoting 
the linear portion of the stress tensor. The quantities c Q jj/p and c a j^p A Bare the 
second- and third-order elastic constants, respectively. Generally, the 
second-order and third-order elastic constants are written in Voigt 
compressed notation whereby a pair of indices is replaced by a single index 
that take on values from 1 to 6 following the notation: 1 1 ->1 , 22^2, 33-»3, 
23->4, 13^5, and 12^6. In equations (5)-(8), T aP , E A b, and w 5 ,y denote the 
biasing stresses, strains, and (static) displacement gradients, respectively. 
Note that the biasing stress T op in the propagating medium is expressed in 
terms of the far-field formation principal stresses (Sv, Sh, and S h ) using 
standard relations that account for stress concentrations caused by the 
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presence of a borehole. In equation (3), AP aj are the perturbations in the 
Piola-Kirchhoff stress tensor elements from the linear portion, P J for the 
reference isotropic medium before the application of any biasing stresses, p Q 
is the mass density of the formation in the reference state, u "} represents the 
eigensolution for the reference isotropic medium for a selected propagating 
mode. The index m refers to a family of normal modes for a borehole in an 
isotropic formation. The quality AP W in equations (5) denotes the increase in 
the wellbore pressure. The frequency perturbations Aw are added to the 
eigenfrequency co m for various values of the wavenumber along the borehole 
axis k z , to obtain the final dispersion curves for the biased state. Note that a 
fractional chance in eigenfrequency from a reference state is equal to a 
fractional change in phase velocity at a given axial wavenumber. Although 
this method of solution is valid for all modes, analyzed herein are stress- 
induced effects only on the Stoneley (m = 0) and flexural (m = 1) modes. The 
volume integrals as shown in equation (3) are computed in the statically 
deformed configuration of the borehole. 

A general perturbation model as defined by equation (3), relates 
perturbations from a chosen reference state in the three principal formation 
stresses (Ao v , Aa H , and Aa h ) to fractional changes in the borehole Stoneley 
and flexural velocities at a given wavenumber kj by 
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where the stress-coefficient of velocity at a given wavenumber kj is given by 
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and the subscript m denotes one of the borehole modes (m=1 -+ fast-flexural; 
m=2 -► slow-flexural; and m=3 -> Stoneley), and the index j represents one of 
the three principal stresses o v , a H , and a h G=1 — ► cj v ; j=2 — ► a H ; and j=3 — ► 
a h ). The superscripts F and S, respectively, denote the fast- and slow- 
flexural wave velocities; and R refers to the flexural wave velocity in the 
chosen reference state. The fast- and slow-flexural waves have radial 
polarization parallel and perpendicular to the maximum horizontal (a H ) stess 
direction, respectively. The four coefficients wj m; w' 2ffy w^. and w^ are 
perturbation integrals evaluated in terms of modal solutions at a given 
wavenumber kj in the reference state. Cm, C155, and c 14 4 are the three 
independent formation nonlinear constants and jj is the shear modulus in the 
reference state. 

The shear modulus in the reference state is determined from the 
Stoneley velocity log calculated at a rather low frequency which is close to the 
tube wave velocity at a reference depth A from the following relationship 



2 



V f,A ' V ST,A 1 

where p f , A is the mud density at depth A; V f , A , and V S T,Aare the mud 
compressional velocity and Stoneley velocity at depth A, respectively. The 
compressional modulus (A + 2u) is obtained from the formation compressional 
velocity at depth A from the relationship 
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A + 2 M = p b , A Vl 4 . (15) 

The dynamic Young's modulus and Poisson's ratio can be readily obtained 
from the Lame parameters (A and n) or the compressional and shear 
velocities in the chosen isotropic reference state. 

Figure 14 is a flow diagram of a routine for practicing an embodiment of 
the invention for determining estimates of maximum horizontal stress 
magnitude S H . This embodiment has particular application when cross-dipole 
waveforms produce wideband dispersions (for example, n the range 1 to 6 
kHz) exhibiting crossovers. A processor, programmed consistent with the 
routine may be, for example, the processor 270 of the surface equipment of 
Figure 3 or a remote processor. If desired, a downhole processor could also 
perform at least part of the technique. Data may, for example, be collected 
and stored using the type of logging apparatus described in conjunction with 
Figures 3-6, although it will be understood that other suitable equipment can 
be utilized. In an example for the present embodiment, a depth interval of 
reasonably uniform lithology is utilized. The borehole diameter, d, is 
measured, for example with a caliper tool. The borehole fluid (mud) 
compressional velocity, V f , is measured, or can be estimated from the mud 
composition, mass density, in-situ pressure and temperature, as is known in 
the art. The formation mass bulk density, p b , and the mud mass density, p f , 
are measured or estimated, as is known in the art. For example, the 
formation mass bulk density can be obtained, for example, from neutron- 
density logging measurements, and the mud mass density can be derived 
using mud weight information from the drilling fluid supplier. 
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The compressional velocity V p and shear velocity V s of the substantially 
undisturbed formation (that is, the relatively far-field region outside any 
mechanically altered annulus) are obtained, for example, from a standard type 
of sonic log (see, for example, Harrison et al., "Acquisition and Analysis of 
Sonic Waveforms From a Borehole Monopole and Dipole Source for the 
Determination of Compressional and Shear Speeds and Their Relation to 
Rock Mechanical Properties and Surface Seismic Data", Society of Petroleum 
Engineers, SPE 20557, 1990). Flexural wave velocities are determined at a 
number of frequencies so that dispersion curves for orthogonal directions can 
be developed at depth levels of interest. 

In the flow diagram of Figure 14, the block 1402 represents selection of 
a depth interval having reasonably uniform lithology, and where cross-dipole 
dispersions exhibit a crossover. The block 1405 represents the measurement 
or estimation of borehole diameter, for example using a caliper tool (in the 
same or a separate logging run). As represented by the blocks 1407 and 
1410 the mud compressional velocity V f , and the mud mass density p f and 
formation mass bulk density p b , are measured or estimated, as previously 
described. The far-field formation compressional velocity V p and shear 
velocity V s are then derived, as represented by the block 1413. These can be 
measured, in conventional fashion, using the sonic logging tool, at a relatively 
long transmitter-to-receiver spacing, to obtain these measurements in the far 
field which are substantially unaffected by damage caused by the drilling of 
the borehole. Alternatively, these and other parameters can be input if they 
are available from another source. The block 1417 represents the 
determination, at the present depth, of the cross-dipole flexural wave 
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velocities at a number of frequencies and the development of the fast and 
slow cross-dipole dispersion curves, in known manner. 

Next, as represented by the block 1420, at the present depth, a 
reference flexural dispersion is computed for an assumed isotropic and 
homogeneous state from the measured V P , monopole shear V s , formation 
bulk mass density p b , mud mass density p f , mud compressional velocity V f , 
and borehole radius a. Reference can be made, for example, to B.K. Sinha, 
"Sensitivity and Inversion of Borehole Flexural Dispersions for Formation 
Parameters", Geophysical Journal international, Vol. 128(1), pp. 84-96, 
January 1997. If the monopole shear velocity is not available from the 
measurements, it can be estimated from the Stoneley (tube-wave) velocity 
(V T ) and mud compressional velocity V f . The reference flexural dispersion is 
for a borehole in an hydrostatically loaded formation subject to a mean 
compressive stress a vo , which is an unknown at this point. Figure 7 illustrates 
decomposition of formation effective stresses in the far-field into a 
hydrostatically loaded reference and perturbations of triaxial stresses. Given 
the overburden stress a v at the chosen depth, a goal of the inversion is to 
estimate the mean hydrostatic stress o vo together with the perturbations in the 
three principal stresses. 

The block 1 425 represents calculating of the sensitivity coefficients 
W\mj, ™imj, w Lj, and w Amj as a function of wavenumber k if for both the fast and 
slow flexural waves in the previously chosen reference state of block 1420. 
These coefficients are calculated in terms of the flexural wave eigensolution in 
the reference state and the near-wellbore stress distributions caused by either 
the maximum horizontal stress parallel (referred to as the fast-flexural) or 
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perpendicular (referred to as the slow-flexural) to the radial polarization of the 
flexural wave. A cost function, e, is defined (block 1430) from differences in 
the fast- and slow-flexural velocities and sensitivity coefficients at various 
wavenumbers k it and known perturbations in stress magnitudes Ao H and Aa h . 
Figure 8 shows the velocity differences as a function of wavenumber that are 
needed in the expression for the cost function z given as: 

1=1 v* 

M 

-Ha -< 3 ) A ^ A -("4 3 -^223)— A °"/, 

-(< 3 -< 3 )— A<7H-{<3-<z)—*<r h > (16) 
M M 

where i=1,2,3, ... n, are the axial wavenumbers (or frequencies) where the 
measured velocity differences (v/ -V/J/V* are calculated. The block 1435 
represents solving for the perturbations in effective stresses Ao H and Aa h that 
would be needed to attain the actual loaded state of the formation in the far- 
field at the chosen depth. In addition to the two stress magnitudes, Cm/p, 
c 155 /p, and C144/M are also solved from a nonlinear least-squares minimization 
of the cost function z as defined by equation (16). Then, as represented by 
the block 1440, the stress coefficients of the fast and slow flexural velocities 
are calculated as a function of axial wavenumber or frequency defined by 

Jf£- = < + < t !f + W ' 3n + w; n ^, (17) 
V; oa v fl fl fi 
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where the superscripts F and S denote the fast- and slow-flexural waves, 
respectively; and the index i refers to the kpth wavenumber. 

Using differences between the measured fast- or slow-flexural and 
reference flexural velocities at various wavenumbers (as shown in Figure 9), 
and the computed sensitivity coefficients (of block 1440) Aov is computed 
(block 1445) from the following relationship: 

The overburden stress S v is calculated (block 1450) by integrating the 
formation mass density from the surface (denoted by 0) to the chosen depth 
(denoted by B), employing the following equation 

S v = f Q p b gdK (24) 
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where p b is the formation bulk density and g is the acceleration due to gravity. 
Then, as represented by block 1455, the effective vertical stress change o v is 
given by 

a v =S v -aP p , (25) 

where a is the Biot parameter for a porous formation, and P P is the pore 
pressure change at the chosen depth B. The mean effective stress o vo 
assumed in the reference state can then be solved (block 1460) using 

a vo =a - Aa v (26) 
Then, estimates of a H and o h are obtained (block 1465) using 

<*a =<7v*+A<7„, (27) 
<T„=<7 vo +Aa A , (28) 

The pore pressures P p at the selected depth are then added (block 1470) to 
obtain 

S v =a v +aP p , (29) 
S H =<j h +aP p , (30) 

S fc =<r fc +oP p , (31) 
The entire procedure can be repeated at other depths (block 1480), and logs 
can be produced. 

Shales typically constitute more than 70% of the formation that a 
wellbore trajectory passes through before reaching the target reservoir. When 
dipole dispersion crossovers are not observed in the overburden shales, 
previously described velocity-gradient and multi-frequency inversion of cross- 
dipole dispersions cannot be used for estimating horizontal stresses. Vertical 
wells in horizontally layered shale, generally, exhibit VTI (transverse isotropy 
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with vertical symmetry axis) anisotropy. Under these circumstances, shear 
isotropy is observed and cross-dipole dispersions nearly coincide with each 
other. This implies that the horizontal stresses are nearly isotropic, at least 
within the radial depth of investigation of the cross-dipole sonic tool and any 
uncertainty in slowness estimates. Figure 10 shows a schematic of such a 
borehole in the presence of an overburden Sv, and equal horizontal stresses 
Sh. However, such shale intervals show compressional velocity and dipole 
gradients with depth. 

Consider a vertical fluid-filled borehole with its axis parallel to the X r 
direction, and isotropic horizontal stresses in the X 2 -X 3 - plane. Biaxial 
formation stresses with a vertical overburden stress as one of the principal 
stresses can be decomposed into a hydrostatically loaded isotropic reference 
and perturbations in the three principal stresses Aov, Aa H , and Aa H as shown 
in Figure 1 1 . Note that the mean stress a V o in the isotropic reference state is 
not known at this point. However, this assumed state is defined from the 
measured compressional and monopole shear velocities at the chosen depth 
so that small perturbations in the three principal stresses Aov, Aqh, and Aqh, 
would lead to the actual in-situ stresses at this depth. All of the stresses in 
Figure 10 are far-field stresses beyond any stress-concentration annulus 
caused by the presence of a borehole. 

A further embodiment of the invention includes determination of the 
horizontal stress S H magnitude in such overburden shales. This technique 
inverts differences in dipole dispersions at the top and bottom of a depth 
interval with a reasonably uniform lithology and mineralogy. 
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Figure 15 is a flow diagram of a routine for practicing the further 
embodiment of the invention. The measurements and other derived inputs 
can be of the same type as previously described. A shale interval with a 
reasonably uniform lithology/mineralogy is identified, as represented by the 
block 1515. 

The block 1520 represents selection of a reference hydrostatically loaded 
state of rock relatively close to the in-situ state of the formation at the top of 
the chosen depth interval (A). The reference state is defined by the 
compressional V P , and shear V s , velocities which are about 1 to 2% less than 
the measured V P and V s at depth A. The formation density p b , mud density 
p f , mud compressional velocity V fj and borehole diameter (2a) are measured 
values at the top of the depth interval A. 

The block 1525 represents calculating of the sensitivity coefficients 
w Lj> w Lj> w Lj>> and W 4mj as a function of wavenumber kj, for both the fast and 
slow flexural waves in the previously chosen reference state of block 1520. 
These coefficients are calculated in terms of the flexural wave eigensolution in 
the reference state and the near-wellbore stress distributions caused by either 
the maximum horizontal stress parallel (again, the fast-flexural) or 
perpendicular (again, the slow-flexural) to the radial polarization of the flexural 
wave. In the case at hand, it is assumed that 0 H = o h . The block 1530 
represents choosing of a general perturbation model relating perturbations 
from a chosen reference state in the three principal formation stresses (Aa v , 
Aa H , and Aa h ) to fractional changes in the borehole flexural velocities 
measured at depth A, and at a given wavenumber kj, by 
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(32) 

where the stress-coefficient of velocity at a given wavenumber kj is given by 

] - W lmi + W L ^<^< ml 03) 



. m i Ji 



and the subscript m denotes one of the borehole modes (m=1 — ► fast-flexural; 
m=2 slow-flexural; and m=3 — ► Stoneley), and the index j represents one of 
the three principal stresses Ao v , Aa H , and Ao h (j=1 — > o v ; j=2 — ► o H ; and j=3 
->a h ). 

Then, as represented by the block 1535, the stress coefficients of the 
flexural velocity as a function of axial wavenumber or frequency are calculated 

as: 

T^=< + <V + ^ +w; "> <34) 
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where the superscripts A denote flexural velocities at the top of the chosen 
depth interval, and the subscript m has been replaced with index i that refers 
to the kj-th wavenumber. A cost function, e, is defined from differences in 
flexural velocities between the bottom (B) and top (A) of the chosen depth 
interval, and sensitivity coefficients at various wavenumbers k } , and unknown 
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perturbations in stress magnitudes Ao H , and Aa h . Figure 12 shows the 
velocity differences as a function of wavenumber that are needed in the 
expression for the cost function e given as: 



-w; n (A<7 v B -A^H^ + wlJlAff" - Act*) 
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M 

where i =1 ,2,3 ... n, are the axial wavenumbers (or frequencies) and where 
the measured velocity differences (v.* -v, A )/v, R are calculated. 

Using a nonlinear least-squares minimization of the cost function, e, the 
block 1545 represents solving for the initial estimates of formation normalized 
nonlinear constants (cm/p, c 155 /p, and cWp) referred to the chosen 
reference state, and the difference in the effective horizontal stresses 
) in this chosen depth interval. Then, using the calculated ratio 

) from block 1545 (or estimated from other sources), the 
block 1550 represents solving the linearized version of equation (37) for the 
unknowns Cm/p, C155/M, and c 14 4/p. Next, the block 1555 represents inverting 
the differences between the measured flexural velocities in the reference state 
and those at depth A as shown in Figure 13, for Acr$ and Ao„ by solving the 
following difference-equations 
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The overburden stress S v can then be calculated (block 1560) by integrating 
the formation mass density from the surface (denoted by 0) to the chosen 
depth (denoted by A) given by the following equation 

Sv = £ P b gdh, (39) 

where p b is the formation bulk density and g is the acceleration due to gravity. 
Then, as represented by the block 1 565, the effective vertical stress change 
Ov at depth A is given by 

a A =S A -aP P A , (40) 
where a is the Biot parameter for a porous formation, and P A is the pressure 
change at the chosen depth A. The mean effective stress o vo assumed in the 
reference state (see Figure 1 1) can then be solved for (block 1570) using 

<7 vo = a A - AaC (41) 

The effective horizontal stress change at depth A, a A , is then determined 
(block 1575) as 

a A = a vo + Aa A (42) 

The pore pressures at the selected depth A can then be added (block 
1 580) to obtain 

S A =(T A +aP A , (43) 
S*=**+aP*, (44) 
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The operations of blocks 1555 through 1580 can then be repeated, as 
represented by the block 1585, to determine an estimate of the horizontal 
stress S B H at depth B by replacing Vj A with Vj B , and S$ with S* . 
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